GIS-CW
Ziyi Chen
2021/1/6
Raw R code and data visualization for CASA0005
The main purpose for this assignment’s report is to find out the spatial distribution of COVID-19 new cases in London borough and to quantify the effect of ‘Lockdown’ policy on the suppress of the increase of new cases based on statistical analysis and spatial autocorrelation.
Data visualization on map
#import library
library(sf)
library(tmap)
library(tmaptools)
library(tidyverse)
library(here)
library(janitor)
library(RColorBrewer)
library(sp)
library(spdep)After completing loading libraries, loading all datasets needed for the analysis:
#read all dataset for the analysis
#London borough shapefile
Londonborough <-
st_read(
here::here(
'data',
'statistical-gis-boundaries-london',
'ESRI',
'London_Borough_Excluding_MHW.shp'
)
) %>%
st_transform(., 27700)## Reading layer `London_Borough_Excluding_MHW' from data source `/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp' using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## projected CRS: OSGB 1936 / British National GridRead csv file–>cumulative new cases in each time period, raw csv files are in Github repo
#cumulative new cases in 3.23-4.22 period
covid_323_422 <-
read_csv(
'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases-323-422.csv',
na = c("NA", "n/a")
) %>%
clean_names()##
## ─ Column specification ─────────────────────────────────────────────────
## cols(
## area_code = col_character(),
## new_cases = col_double(),
## `population(10k)` = col_double(),
## `new_cases per 10k population` = col_double()
## )#cumulative new cases in 5.23-6.22 period
covid_523_622 <-
read_csv(
'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases-523-622.csv',
na = c("NA", "n/a")
) %>%
clean_names()##
## ─ Column specification ─────────────────────────────────────────────────
## cols(
## area_code = col_character(),
## new_cases = col_double(),
## population_10k = col_double(),
## `new_cases_per_10k_523-622` = col_double()
## )names(covid_323_422)## [1] "area_code" "new_cases"
## [3] "population_10k" "new_cases_per_10k_population"names(covid_523_622)## [1] "area_code" "new_cases"
## [3] "population_10k" "new_cases_per_10k_523_622"Connect the CSV file data and shapefile
#join data with shapefile and new cases data in 3.23-4.22 period
Londonborough_merged <-
Londonborough %>% left_join(covid_323_422, by = c('GSS_CODE' = 'area_code')) %>%
distinct(GSS_CODE, NAME, new_cases,new_cases_per_10k_population) #choose which column to be used in the following analysis
#join data with shapefile and new cases data in 5.23-6.22 period
Londonborough_merged_523 <-
Londonborough %>% left_join(covid_523_622, by = c('GSS_CODE' = 'area_code')) %>%
distinct(GSS_CODE, NAME, new_cases, new_cases_per_10k_523_622) #choose which column to be used in the following analysis#making maps
tmap_mode('view')## tmap mode set to interactive viewingbreaks_323 = c(10, 15, 20, 25, 30, 35)
tm1 <- tm_shape(Londonborough_merged) +
tm_polygons(
'new_cases_per_10k_population',
breaks = breaks_323,
alpha = 0.5,
palette = 'PuBu',
colorNA = 'white',
title = 'New cases rate'
) + tm_layout(
legend.position = c('left', 'bottom'),
legend.outside = FALSE,
legend.title.size = 1.2,
legend.text.size = 0.75,
frame = FALSE
) + tm_credits('(A) New cases rate in 3.23-4.22', position = c('left', 'top'), size = 1.2) + tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =c('right', 'bottom'), text.size = 0.6)
#view map 1
tm1## Credits not supported in view mode.## Compass not supported in view mode.## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.#saving map
tmap_save(tm1,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(323-422).png')## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(323-422).png## Resolution: 2389.896 by 1845.268 pixels## Size: 7.966321 by 6.150895 inches (300 dpi)#making map 2
breaks_523=c(0,0.5,1,1.5,2,2.5,3)
tm2 <- tm_shape(Londonborough_merged_523) +
tm_polygons(
'new_cases_per_10k_523_622',
breaks = breaks_523,
alpha = 0.5,
palette = 'PuBu',
colorNA = 'white',
title = 'New cases rate'
) + tm_layout(
legend.position = c('left', 'bottom'),
legend.outside = FALSE,
legend.title.size = 1.2,
legend.text.size = 0.75,
frame = FALSE
) +
tm_compass(type = "arrow", position = c("right", "bottom")) + tm_credits('(B) New cases rate in 5.23-6.22',
position = c('left', 'top'),
size =
1.2) +
tm_scale_bar(position = c('right', 'bottom'), text.size = 0.6)
tm2## Credits not supported in view mode.## Compass not supported in view mode.## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.tmap_save(tm2,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(523-622).png')## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(523-622).png## Resolution: 2389.896 by 1845.268 pixels## Size: 7.966321 by 6.150895 inches (300 dpi)#combine two map
t=tmap_arrange(tm1,tm2,nrow=1)
tmap_save(t,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/New_cases_rate.png')## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/New_cases_rate.png## Resolution: 2100 by 2100 pixels## Size: 7 by 7 inches (300 dpi)t## Credits not supported in view mode.## Compass not supported in view mode.## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.## Credits not supported in view mode.## Compass not supported in view mode.## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.###Global Moran
#data cleaning for 3.23-4.22 data
london_exclude_city <-
na.omit(Londonborough_merged) #exclude rows contain null data
#data cleaning for 5.23-6.22 data
london_exclude_city_523 <- na.omit(Londonborough_merged_523)
#creating polygon for 323 dataset and 523 dataset
neibour_323 <- poly2nb(london_exclude_city, queen = TRUE)
neibour_323[[1]]## [1] 19 20 21 22neibour_523 <- poly2nb(london_exclude_city_523,queen=TRUE)
neibour_523[[1]]## [1] 19 20 21 22#assign weight matrix for each neighbouring polygon using row standardisation
lw_323 <- nb2listw(neibour_323, style = "W", zero.policy = TRUE) #each neighbour is assigned a quarter of total weight
lw_523 <- nb2listw(neibour_523,style='W',zero.policy = TRUE)
#computing global Moran's I
moran.test(london_exclude_city$new_cases_per_10k_population, lw_323)##
## Moran I test under randomisation
##
## data: london_exclude_city$new_cases_per_10k_population
## weights: lw_323
##
## Moran I statistic standard deviate = 3.1736, p-value = 0.0007528
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.35308476 -0.03225806 0.01474307moran.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)##
## Moran I test under randomisation
##
## data: london_exclude_city_523$new_cases_per_10k_523_622
## weights: lw_523
##
## Moran I statistic standard deviate = 2.5898, p-value = 0.004802
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.28295263 -0.03225806 0.01481388###Local Moran visualization
#local moran value for data 3.23-4.22
local_moran_323 <- localmoran(london_exclude_city$new_cases_per_10k_population, lw_323)
summary(local_moran_323)## Ii E.Ii Var.Ii Z.Ii
## Min. :-0.419544 Min. :-0.03226 Min. :0.1102 Min. :-0.71017
## 1st Qu.:-0.002326 1st Qu.:-0.03226 1st Qu.:0.1664 1st Qu.: 0.08039
## Median : 0.149082 Median :-0.03226 Median :0.2155 Median : 0.36073
## Mean : 0.353085 Mean :-0.03226 Mean :0.2460 Mean : 0.82412
## 3rd Qu.: 0.619470 3rd Qu.:-0.03226 3rd Qu.:0.2974 3rd Qu.: 1.59790
## Max. : 2.209853 Max. :-0.03226 Max. :0.4612 Max. : 4.11139
## Pr(z > 0)
## Min. :0.0000197
## 1st Qu.:0.0550665
## Median :0.3595183
## Mean :0.3109317
## 3rd Qu.:0.4679710
## Max. :0.7612010#local moran value for data 5.23-6.22
local_moran_523 <- localmoran(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
summary(local_moran_523)## Ii E.Ii Var.Ii Z.Ii
## Min. :-0.40785 Min. :-0.03226 Min. :0.1106 Min. :-0.80748
## 1st Qu.:-0.05669 1st Qu.:-0.03226 1st Qu.:0.1670 1st Qu.:-0.05456
## Median : 0.16490 Median :-0.03226 Median :0.2164 Median : 0.43173
## Mean : 0.28295 Mean :-0.03226 Mean :0.2470 Mean : 0.68948
## 3rd Qu.: 0.44756 3rd Qu.:-0.03226 3rd Qu.:0.2986 3rd Qu.: 1.17418
## Max. : 1.39046 Max. :-0.03226 Max. :0.4632 Max. : 3.48158
## Pr(z > 0)
## Min. :0.0002492
## 1st Qu.:0.1213157
## Median :0.3343064
## Mean :0.3309740
## 3rd Qu.:0.5217552
## Max. :0.7903052#plot local moran map
#There are 5 columns of data.
#copy some of the columns (the I score (column 1) and the z-score standard deviation (column 4))
#the z-score (standardised value relating to whether high values or low values are clustering together)
#change local_moran type
local_moran_tibble_323 <- as_tibble(local_moran_323) #change to dataframe
local_moran_tibble_523 <- as_tibble(local_moran_523)
london_exclude_city <-
london_exclude_city %>% mutate(local_moran_I = as.numeric(local_moran_tibble_323$Ii)) %>% mutate(local_moran_z =as.numeric(local_moran_tibble_323$Z.Ii))
london_exclude_city_523 <-
london_exclude_city_523 %>% mutate(local_moran_I = as.numeric(local_moran_tibble_523$Ii)) %>% mutate(local_moran_z =
as.numeric(local_moran_tibble_523$Z.Ii))
#setting color
MoranColours <- rev(brewer.pal(8, "RdBu"))
#plot a map (in Rmd documend, tmap mode changes to viewing, but raw code is plotting mode)
tmap_mode('view')## tmap mode set to interactive viewing#plotting local moran map for 3.23-4.22
tm_323 <- tm_shape(london_exclude_city) +
tm_polygons(
'local_moran_I',
alpha = 0.5,
palette = MoranColours,
title = 'Local Moran I',midpoint=NA
) + tm_layout(
legend.position = c('left', 'bottom'),
legend.outside = FALSE,
legend.title.size = 1.2,
legend.text.size = 0.75,
frame = FALSE
) + tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =
c('right', 'bottom'), text.size = 0.6)+tm_credits('(A) Local Moran in 3.23-4.22',position = c('left','top'),size = 1.2)
tm_323## Credits not supported in view mode.## Compass not supported in view mode.## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.tmap_save(tm_323,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_323.png')## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_323.png## Resolution: 2389.896 by 1845.268 pixels## Size: 7.966321 by 6.150895 inches (300 dpi)#plotting local moran map for 5.23-6.22 dataset
tm_523 <- tm_shape(london_exclude_city_523) +
tm_polygons(
'local_moran_I',
alpha = 0.5,
palette = MoranColours,
title = 'Local Moran I',midpoint=NA
) + tm_layout(
legend.position = c('left', 'bottom'),
legend.outside = FALSE,
legend.title.size = 1.2,
legend.text.size = 0.75,
frame = FALSE
) +tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =
c('right', 'bottom'), text.size = 0.6)+tm_credits('(B) Local Moran in 5.23-6.22',position = c('left','top'),size = 1.2)
tm_523## Credits not supported in view mode.## Compass not supported in view mode.## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.tmap_save(tm_523,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_523.png')## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_523.png## Resolution: 2389.896 by 1845.268 pixels## Size: 7.966321 by 6.150895 inches (300 dpi)#combine two map
t_c=tmap_arrange(tm_323,tm_523,nrow=1)
t_c## Credits not supported in view mode.## Compass not supported in view mode.## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.## Credits not supported in view mode.## Compass not supported in view mode.## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.tmap_save(t_c,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_combined.png')## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_combined.png## Resolution: 2100 by 2100 pixels## Size: 7 by 7 inches (300 dpi)###Geary’s test result
#geary's test
#geary's test for 3.23-4.22 data
geary.test(london_exclude_city$new_cases_per_10k_population, lw_323)##
## Geary C test under randomisation
##
## data: london_exclude_city$new_cases_per_10k_population
## weights: lw_323
##
## Geary C statistic standard deviate = 2.9614, p-value = 0.001531
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.63209605 1.00000000 0.01543425#geary's test for 5.23-6.22 data
geary.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)##
## Geary C test under randomisation
##
## data: london_exclude_city_523$new_cases_per_10k_523_622
## weights: lw_523
##
## Geary C statistic standard deviate = 2.1708, p-value = 0.01497
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.73102781 1.00000000 0.01535192###Getis Ord General G result
#Getis Ord General G for 3.23-4.22
globalG.test(london_exclude_city$new_cases_per_10k_population,lw_323)## Warning in globalG.test(london_exclude_city$new_cases_per_10k_population, :
## Binary weights recommended (sepecially for distance bands)##
## Getis-Ord global G statistic
##
## data: london_exclude_city$new_cases_per_10k_population
## weights: lw_323
##
## standard deviate = 2.5827, p-value = 0.004901
## alternative hypothesis: greater
## sample estimates:
## Global G statistic Expectation Variance
## 3.347784e-02 3.225806e-02 2.230541e-07#Getis Ord General G for 5.23-6.22
globalG.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)## Warning in globalG.test(london_exclude_city_523$new_cases_per_10k_523_622, :
## Binary weights recommended (sepecially for distance bands)##
## Getis-Ord global G statistic
##
## data: london_exclude_city_523$new_cases_per_10k_523_622
## weights: lw_523
##
## standard deviate = 1.525, p-value = 0.06362
## alternative hypothesis: greater
## sample estimates:
## Global G statistic Expectation Variance
## 3.336104e-02 3.225806e-02 5.230748e-07###Descriptive statistics summary
library(ggplot2)
#histrogram
#only use rate of new cases to plot histrogram
covid_323_hist <- hist(covid_323_422$new_cases_per_10k_population,col='Dark blue',density=10,xlab='Rate of new cases',main='Rate of new cases during 3.23-4.22',angle=45,ylim=c(0,15))
mean_323 <- mean(covid_323_422$new_cases_per_10k_population)
var_323 <- var(covid_323_422$new_cases_per_10k_population)
std_323 <- sd(covid_323_422$new_cases_per_10k_population)
mean_323## [1] 23.5875var_323## [1] 25.43532std_323## [1] 5.043344covid_523_hist <- hist(covid_523_622$new_cases_per_10k_523_622,col='Dark blue',density=10,xlab='Rate of new cases',main='Rate of new cases during 5.23-6.22',angle=45,ylim=c(0,15))
mean_523 <- mean(covid_523_622$new_cases_per_10k_523_622)
var_523 <- var(covid_523_622$new_cases_per_10k_523_622)
std_523 <- sd(covid_523_622$new_cases_per_10k_523_622)
mean_523## [1] 1.771875var_523## [1] 0.2949899std_523## [1] 0.5431297